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Abstract. We study large fluctuations in evolutionary games belonging to the coordi- 
nation and anti-coordination classes. The dynamics of these games, modeling cooperation 
dilemmas, is characterized by a coexistence fixed point separating two absorbing states. We 
are particularly interested in the problem of fixation that refers to the possibility that a 
few mutants take over the entire population. Here, the fixation phenomenon is induced by 
large fluctuations and is investigated by a semi-classical WKB (Wentzel-Kramers-Brillouin) 
theory generalized to treat stochastic systems possessing multiple absorbing states. Impor- 
tantly, this method allows us to analyze the combined influence of selection and random 
fluctuations on the evolutionary dynamics beyond the weak selection limit often considered 
2j^' in previous works. We accurately compute, including pre-exponential factors, the probability 

distribution function in the long-lived coexistence state and the mean fixation time necessary 
for a few mutants to take over the entire population in anti-coordination games, and also 
\^ ' the fixation probability in the coordination class. Our analytical results compare excellently 

fSJ , with extensive numerical simulations. Furthermore, we demonstrate that our treatment is 

superior to the Fokker-Planck approximation when the selection intensity is finite. 
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I. INTRODUCTION & MODELS 

Systems in which successful strategies spread by imitation or reproduction can be described by evolutionary game 
theory (EGT), whose prototypical models are commonly studied in the context of evolutionary biology, ecology, 
sociology and economics [2-||. Recently, it has also been realized that techniques of statistical physics can help 
gain further insight into this interdisciplinary area [0, R] . Originally, EGT was formulated in terms of deterministic 
replicator equations valid to treat populations of infinite size [^ g, || and related to the celebrated Lotka-Volterra 
equations |j5-0, |rpO|- However, it has been long recognized that the picture emerging from replicator dynamics is 
often fundamentally altered by demographic fluctuations. 



One of the most striking effects of fluctuations in EGT is fixation that refers to the possibihty for a few "mutants" 
to take over (fixate) an entire population causing the extinction of the "wild species" . To rationalize the effect of 
stochasticity in finite populations, Nowak et al. [g, Q introduced a parameter w controlling the interplay between 
random (demographic) fluctuations and selection (leading to nonlinear effects). This approach reflects the fact that 
evolutionary processes comprise two main competing mechanisms. On the one hand, there is selection by individuals' 
different fitness (reproductive potential), which underlies adaptation [0, || |5| ^, jf^-jl^l- On the other hand, in all 
finite-size populations birth and death events cause random (demographic) fluctuations, which play a central role in 
neutral theories where selection is considered of marginal importance [|l6|-|l9| . 

Most of the analytical results concerning EGT in finite population have been obtained in the weak selection limit 
{w — >■ 0), which is often biologically relevant |2^ and greatly simplifies the mathematical analysis. In particular, 
the fixation probability of a species under frequency-dependent selection in a finite population (of two species) has 
been computed in this limit |ll|. However, very different behaviors have been obtained under strong and weak 
selection (see, e.g., pl|-25t) and the respective influence of selection and demographic fluctuations on evolution still 
remains to be fully understood. Our purpose in this paper is to study fixation resulting from large fluctuations in 
two classes of evolutionary games (the anti-coordination and coordination classes, see below) under arbitrary (finite) 
selection strength, and to elucidate the nontrivial interplay between selection intensity and demographic fiuctuations. 

For the sake concreteness, throughout this paper we will consider (symmetric) 2x2 evolutionary games. Here, one 
has a homogeneous (well-mixed) population comprising a total of N individuals, n are of the "mutant" species A and 
N — n are of the "wild" type B. As usual, it is assumed that there are pairwise and symmetric interactions between 
individuals drawn at random. The reproductive potential of an individual is specified by the payoff of interaction with 
others |Q-|[ |[. Specifically, when two individuals of species A interact, both receive a payoff a. The interaction of a 
pair of individuals of different species yields payoffs b and c to the A and B individuals, respectively. Similarly, when 
two individuals of type B interact, both get a payoff d. Therefore, the respective average (per individual) payoffs of 
species A and B read: HA(n) = {n-l)/{N-l)a + {N-n)/{N -l)b, andHij(n) = n/{N-l)c+{N-n-l)/{N-l)d, 
with self-interactions being excluded. It is useful to introduce the difference between the average payoffs AIl{n) — 
H^(n) - Usin) = {a - b ~ c + d)n/{N - 1) - {d - b)N/{N - 1) - (a - d)/{N - 1). For populations of infinite size 
{N — )> oo), the dynamics is of mean-field type and commonly described by the replicator equations. The latter are 
obtained from the payoff matrix by comparing the success of a given type with the population average flj-l^l ■ In this 
limit X = n/N can be regarded as a continuous variable and the dynamics is specified by the following replicator 
equation: 

X = x{l - x)An{x) , An{x) ^{a-c)x-{d-b){l-x). (1) 

This equation is characterized by two absorbing fixed points x = x\ = 1 and x = x*^ = corresponding to a system 
with all A's and B's, respectively. Moreover, there can be an interior fixed point obtained by solving AH(x*) — 

'^-^ (2) 



a — b — c + d 

Depending on the entries of the payoff matrix, one has various classes of games representing models of cooperation 
dilemmas [0, O, E5J. When one species always dominates the other, one has the dominance class, where A is the 
dominant type when a > c and b > d, while B is dominant when a < c and b < d. In this work we are interested in 
the other two classes, anti-coordination games (ACG) and coordination games (CG), where there is an interior fixed 
point corresponding to a coexistence state. In the class of ACG {e.g., "snowdrift" and "hawk-dove" games |l], |[ |[), 
c > a, b > d, and x* is an attractor corresponding to the stable coexistence of A and B types. Here, the absorbing 
states cc^ and x*^ are (evolutionary) unstable. On the other hand, in the class of CG {e.g. "stag-hunt" game |l]-||) 
a > c, d > b, and x* is repelling, while the absorbing states are (evolutionary) stable. 

Fluctuations arising from the discreteness of individuals and from the stochastic nature of the interactions may 
drastically alter the predictions of the deterministic replicator equations |26|| . In particular, a few mutants can always 
attain fixation by taking over the entire population (see below). The resulting stochastic evolutionary dynamics is 
aptly described in terms of single-step birth-death processes P7|-p9[ , for which a key quantity is the time-dependent 
probability distribution function (PDF) P„(i) of population sizes. The latter gives the probability of finding the 
system in a state with n individuals of species A at time t, and obeys the following master equation: 

^^^ = T+{n - l)P„_i +T-{n + l)P„+i - [T+{n) + T-{n)]P,,. (3) 

Here T+(n) and T^{n) respectively denote the transition rates from a state with n A's to a state with n + 1 and 
n— \ A's. As the state space is bounded, n E [0, N], and n = and n = N are absorbing states, the transition rates 
at the boundaries satisfy T±(0) = T^{N) = 0. 



According to general prescriptions of EGT, the transition rates are functions of each species' fitness (effective 
potential to reproduce), f„, with a G (A, B), i.e. T'^{n) = T'^[f„{n)]. The fitness of an individual of species a 
reads ||, f^{n) = I — w + wlla{n). Here, the interplay between random fluctuations and selection is controlled 
by the parameter w (with < w < 1), where in the neutral case, w = 0, there is no selection (but only random 
fluctuations), while in the strong selection regime, w = 1, the influence of random fluctuations is negligible. 

In this paper we consider the following update rules commonly used in the EGT literature, specifying the functional 
dependence of T^ on the species fitnesses. For the frequency-dependent Moran Process (fMP) [0, 0, Bfl, one has 






T {n) = -=— -$(n). 



(4) 



where f{n) = [nfA{n) + (N — n)fB{n)]/N is the population average fitness and $(n) = n{N — n)/N'^. For the linear 
Moran process (LMP) [Q ^, the transition rates read 



T+in) = i{l + IfAin) - fin)]} ^n) , T-{n) = ^{1 + [/^(n) - ./»]}$(«). 



(5) 



As the LMP is obtained from a small w expansion of the fMP, (||) can be regarded as the "weak selection" counterpart 
of the rates (^. A process closely related to the LMP is the "local update" process (LUP) [|o| ^ with rates: 

T+{n) = ^{l + [fAin)~fB{n)]}Hn), T-(n) - ^{1 + [/^(n) - /^(n)]} <l>(n). (6) 

Here, for simplicity and without restriction, we have assumed that the maximum payoff difference is 1 pOj . Finally, 
for the Fermi process (FP) [|, H H, |§, one has 

r+(n) = {I + exp[/B(n)-/^(n)]}-^$(n), r-(n) = {1 + cxp[/^(n) - /^(n)]}-^ $(n). (7) 



In the following, we omit in all these cases the self- interaction terms in the expressions of Ha(".) and Hs(n) [ p3[ . 
Note, that multiplying both sides of (0) by n and summing over all n's, one obtains an equation for the mean number 
of A individuals. In the leading order of A^ S> 1 and upon rescaling time, one arrives at the following rate equation 
for the concentration of mutants: x = T^[Nx) — T^ (Nx). Such a replicator-like equation shares the same properties 
(fixed points and stability) of Eq. (|^), but generally differs from it when / is non-constant (see e.g., [^ p2|). 

Let us denote by (j)f the probability that i mutants of species A (usually i <^ N) replace all the individuals of the 
wild type B. That is, (^^ is the probability of fixation of the A species starting with i mutants. The conditional and 
unconditional mean fixation times (MFTs) t^ and r^, respectively, are the times associated with the fixation event. 
The former, rf', is the average time it takes for i mutants of species A to take over the population, while the latter, r^, 
is the mean time it takes the population, initially comprising i individuals of species A, to become homogeneous again 
(ie.populated either by all A's or all B's). For all one-dimensional single-step birth-death processes, as those defined 
by (^0), there are exact formulas for the above quantities |g^, ^ |7| ^. For instance, the fixation probability reads 



i+E^\nti 
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where 7i = T^{i)/T^{i), while for the unconditional MET, one has 



(8) 
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(9) 



Even though the expressions d^^) are exact, they are unwieldy and cannot be generalized to multi-step processes and 
to s X s games, with s > 2. Furthermore, it is highly nontrivial to extract their asymptotic behavior. In fact, with 
the exception of Ref. ||2J| where the fixation probability and MFTs were calculated in the leading order for the fMP 
(focusing on w — 1), most of the analytical results in the literature have been obtained in the weak selection limit, 
Nw < 1, often using the Fokker-Planck equation (EPE) HJ 



25, 30, 31 



In this paper we go beyond the weak selection limit and investigate fixation phenomena induced by large fluctuations 
in the classes of ACQ and CG. Our approach relies on the WKB approximation |3J] applied directly to the master 
equation (0) [p5|-p8[, that we generalize to account for the existence of multiple absorbing states. Here, the WKB 
approximation is an asymptotic series expansion in powers of 1/A'^ <C 1 based on an exponential ansatz |39 [see 
Eq. (O)] and differs from the van-Kampen system size expansion that yields the FPE p^, S (see Sec. III.B). With 



the WKB approach and for any finite selection intensity, we accurately compute the MFTs and fixation probability 
(including pre-exponential factors) for generic transition rates. Our general results are then applied to the processes 
defined by the transition rates (Qj^) ^^'^ successfully compared with extensive numerical simulations. The predictions 
of our WKB approach are also shown to be advantageous over those of the FPE when w is finite (see also Refs. pof). 
The remainder of this paper, of which a brief account has recently been given in Ref. M], is organized as follows. 
The next section is dedicated to the ACG, for which a general treatment is presented in Section II. A, while applications 
are discussed in Section II. B. Section III concerns the CG class, with general results and applications discussed in 
Sections III. A and III.B, respectively. Finally, we give a summary of our findings and present our conclusions in 
Section IV. Some technical details are relegated in an appendix. 

II. ANTI-COORDINATION GAMES: METASTABILITY, FIXATION TIMES AND PROBABILITY 

In this section wc use the WKB approach to investigate large fluctuations in systems of evolutionary games char- 
acterized by metastability. For iV ^ 1, in general this case arises in the ACG, to which the snowdrift and hawk-dove 
games belong [y-Q. In ACG, the elements of the payoff matrix satisfy b > d and c > a. Here, the attracting (interior) 
fixed point x* [see Eq. (||)] in the language of the replicator equation (|l|), separates the repelling fixed points x\ (all 
A's) and Xg (all B's). In the presence of internal noise, x* corresponds to a long-lived metastable state, where after a 
sufficiently long time the system is eventually driven into one of the two absorbing states via a large fluctuation. In 
this section we first derive general results concerning the metastable dynamics of stochastic systems possessing two 
absorbing states. Then, using the transition rates (H)-(n)i '^s apply these findings to the case of ACG. 

A. General treatment & results 

Our starting point is the master equation (||). In case of a finite space n G [0, A^], one can always expand P„(t) in 
a finite series of eigenvectors and eigenvalues of the stochastic generator associated with the Markov chain (^) . We 
assume here and henceforth that A'^ ^ 1. In this case, for any (sufficiently large) given initial population of A's, after 
a relaxation time scale t^, the system converges into a long-lived metastable (coexistence) state prior to fixation of 
either species. This metastable state corresponds to a PDF of population sizes peaked in the vicinity of n* = Nx*, 
where x* is attracting interior fixed point of (|l|). Here, the MFT r is associated with the slow decay of the metastable 
state, see below, and satisfies U <C r. That is, fixation occurs in the aftermath of a long-lasting coexistence state. 

It turns out, that for times t ^ t^, when the system has already converged into the metastable state, the higher 
excited eigenvectors have already decayed (see, e.g., Ref. |Q). At such times, only the first excited eigenvector of 
(0), 7r„, called the quasi-stationary distribution (QSD), survives and determines the shape of the metastable PDF. 
Correspondingly, the decay rate of the metastable PDF is determined by the first-excited eigenvalue of (^ [^ . While 
the metastable PDF decays, the probabilities Po{t) and PN{t) slowly grow in time and, at t — > oo, reach some final 
values such that Po(oo) -I- Pi\[{oo) = 1. Therefore, at i ^ tr, one can write 

Pn=i,...M~i{t)^7Tr.e''/^ , Po (*) ^ <^(1 " e"*/^) , P^(i) ~ (1 - 0) (1 - e"*/^). (10) 

From this metastable dynamics one can immediately see that r is the (unconditional) MFT, while cj) — (p^ is the 
fixation probability of species B, and 1 — is the fixation probability of species A. 

It follows from (nO|) that Pq + Pjv — (l/''')e~*^^j while from the master equation (^ one has Pjv = T^{N — 
l)7rAr_ie-*/^ and Pq = p- {l)nie^*/^ . We thus obtain Po{t) = [r^i - r+(A^ - l)7rAr_i] e"*/^ = T-{l)nie-*/^ , 
whose solution (with Po(0) = 0) is Po{t) = [I - tT+{N - l)7rAr_i] (1 - e"*/"'). Using this solution and Eq. (|l^), we 
obtain the fixation probability which turns out to be the relative flux into the absorbing state n = 0. Moreover, as 
the unconditional MFT is the (inverse of the) decay rate of the metastable state, it is given by the (inverse of the) 
sum of probability fluxes into the two absorbing states . Thus, we have 

= r-(l)^ir, r=[p-(l)^i+r+(A^-l)^jv-i]-\ (11) 

where the unknowns tti and Tr^r^i will be determined shortly. Correspondingly, r and r^ - the conditional MFTs 
of species A and B respectively, can also be found. The former is the mean time it takes the A species to flxate 
conditioned on the non-flxation of the B species; it is determined by the the inverse of the flux to the state n = N. 
Using the same reasoning for t^ , we thus have r"^ = [T+(A^ — l)7rjv-i]^^ and t^ = [P^(1)7Ti1^^. Note, that since 
we have assumed that the system reaches the metastable state prior to fixation, our results (O) are independent of 
the initial condition. As we shall see below, this assumption holds for ACG when the selection strength is finite. 



Substituting the metastable ansatz ( JlOl ) into Eq. (|3|), we arrive at the quasi-stationary master equation (QSME). 
Neglecting the exponentially small term 7r„/T (to be verified a-posteriori) on the left-hand-side, the QSME becomes 



^ T+{n - l)7r„_i + r-(n + l)7r„+i - [T+{n) + T-{n)]^„ 



(12) 



In the following, this equation is analyzed by using the WKB approximation. Our aim, in addition to finding the 
fixation probability and MFTs, is to find the complete QSD, 7r„, and demonstrate the non-Gaussian nature of its 
tails. To do so, and since there are two absorbing states in this problem, we solve the QSME (|l3) separately in three 
overlapping regions: in the bulk and not too close to the absorbing boundaries, and in the close vicinities of n = 
and n ~ N . These solutions are then matched in their regions of joint validity. 

In the bulk region (whose accurate boundaries are specified below) we employ the WKB ansatz 



^,^ = n{x) = ^e-^^(-)-^^(-) , 



(13) 



where A^ ^ 1, and we have introduced the rescaled coordinate x ^ n/N. Here, S{x) is the action while Si{x) is the 
amplitude. To have a consistent perturbation theory, we assume that these quantities are smooth functions of order 
unity. The constant prefactor A is introduced for technical convenience. It is convenient to rewrite the transition 
rates as continuous functions T±{x) = T^{n) of the rescaled continuous coordinate < a; < 1. We will assume that 
in the bulk, T±{x) — 0(1), which is satisfied by all the transition rates (0)-(n). 

Plugging ansatz (O) into Eq. ([l2|), and expanding the functions of x ± N^-^ up to 0{N^^), we arrive at 



7r(a;)<^r+(.T) 



.s' 



1-^ + ^1-1 

2N N 



T-ix) 



^ 2N N ^ 



1 

N 



^-s' 



TL{x)~e'Ti{x) =0-(14) 



This equation can be solved order by order in N '^ 1. In the leading 0(1) order, one obtains a stationary Hamilton- 
Jacobi equation for the action S{x), H[x, S'{x)] — with a Hamiltonian given by 



Hix.p) = T+{x)ieP - 1) + T-ix)ie-P - 1) . 



(15) 



where we have defined the auxiliary momentum coordinate p{x) = S'{x) |36|] . Therefore, the leading-order calculations 
correspond to finding a nontrivial zero-energy trajectory of the auxiliary Hamiltonian (051) ISO]. It turns out that 
there is exactly one (real) such trajectory [|8J, Paix), called the activation trajectory. It represents the "optimal-path" 
along which the stochastic system evolves, almost with certainty, from the metastable state towards fixation. Here, 
the solution of H[x,Paix)\ = is Pa{x) = — \n[T+{x)/T-{x)] [pSi. The corresponding action along this trajectory is 



S{x) 



HT+io/T-mdc- 



(16) 



In the subleading 0(1/A^) order, one arrives at a first-order transport-like differential equation for Si{x), whose 
solution is [^ ISSl 



Si{x)^^HT+{x)T-{x)] 



Therefore, the solution in the bulk region can be written as 



7r(a;) 



A 



^n{x)T-{x) 



^Nr\n[r+iO/T-(.0]di 



(17) 



(18) 



It is worth emphasizing that this solution is valid in the regime where 7±(x) = 0(1), i.e., not too close to the absorbing 
boundaries a; = and x — 1, where the transition rates vanish. As for A'^ ^ 1 the QSD is strongly peaked in the vicinity 
of the attracting fixed point x* , the constant A can be determined by normalizing to unity the Gaussian asymptote 
of the QSD around x*. The latter is obtained by expanding Eq. (Ol) to second order about x* and using Paix*) = 
[since T+{x*) — T-{x*)]. Normalizing the resulting Gaussian asymptote, 7r(a;) ~ j[g-NSix*)-Si{x*)-{N/2)s"ix-'){x-x') ^ 
yields the constant A, and therefore the QSD is given by 



7r(x) = T+{x*) 



S"{x*) 



-N[S(x)-S{x*)] 



(19) 



^ 2T:NT+ix)T-{x) 
Here, S"{x*) = TL{x*)IT-{x*) - T[{x*)/T+{x*) > 0, as x* is an attracting fixed point and so Tlix*) - TLix*) < 0. 



We now turn to dealing with the QSME (|lj) in the close vicinities of the absorbing boundaries where the WKB 
approximation breaks down. First, expanding the transition rates T±{x) ~ xT^{Q) about x — Q [where 7±(0) ~ 0], 
and multiplying Eq. (P) by N, one has = T|(0)(n - l)7r„_i + rZ(0)(n + l)7r„+i - n[T+{Q) + TZ(0)]7r„, whose 
recursive solution is Pq] 

^. = ^-P^-^. (20) 

(i?o-l)n ^ ' 

Here we have introduced the parameter Rq = T^{0)/TL{0)- This procedure turns out to be valid in the range 



1 < n <^ Vn IQ. Similarly, close to the boundary n — N, the rates in the QSME (12) can be expanded in the 
vicinity of x = 1 yielding = Tlil){N - n + l)7r„_i + T1{1){N -n~ l)7r„+i - {N ~ n)[r|(l) + TLil)]-Kn- The 
solution of this equation, valid for 1 < iV — rt <C \N and with i?i = 7^(1)/7^(1), satisfies 



{Ri~l){N^n) 



We are in the position to find the complete QSD by matching Eq. (pQ) and (p^) with the asymptotes of (nS) in the 
regions 1 ^ n <C VN and 1 <^ N — n <^ Vn, respectively. In the vicinity of cc = 0, the asymptote of ([19^ can be 
found by writing S{x) ~ 5(0) + xpa{0), with pa{0) = ln[rZ(0)/r+(0)], which yields 

Trfx) = '^+(^*)V^'^"(^*) j^Nx^-N[SiQ)~Six-')] /22) 

^ ' x^27rNTliO)Tl{0) " ■ ^ ' 



This asymptote is valid for 1 ^ n ^ v A^ ||38| . Similarly, the asymptote of Eq. ( |19| ) in the vicinity of a; = 1 reads 

^ r+(x*) vy'(x*)<<^'"^ g-^[s(i)-s(.*)] (23) 

and is valid for K A^ - n < %/iV. Matching Eqs. (|2|) and (||), respectively with Eqs. (|^ and ^ yields 



V 2^ ^/UWiW) ' V 27r ^r;(i)rz(i) ' ^ ' 

With the expressions ( |24| ) and (|9|), the QSD has been completely determined. The fixation probability and the 
MFTs can then be computed according to (|ll|). In fact, as r_(l) ~ {l/N)TL{Q) and T+{N - 1) ~ (l/7V)|r+(l)| [as 
7^(1) < 0], the fixation probability and unconditional MET (pl| ) become 

-^'(^^-^ - ^ (25) 



" r:(o)7ri + |r;(i)|^jv-i ' TL{Q)ni + |r;(i)|^iv-i ' 

while the conditional MFTs of species A and B are respectively t^ = N[\T+{l)\iTN~ir^ and r^ = iV[Tl(0)7ri]"^ 
Importantly, since we have assumed that r is exponentially large in N , these results indicate that our theory is valid 
as long as N[S{1) - S{x*)] > 1 and A^[5(0) - S{x*)] > 1. 

B. Applications 

As applications of the general results that we have derived, we now explicitly consider the four types of update 
rules mentioned above, i.e., the fMP, LMP, LUP and the FP (|)-(0). For each of them we obtain the QSD, the MFTs 
and the fixation probability under arbitrary (but finite) selection strength < w < 1. 

1. Frequency-dependent Moran process 

For the fMP, the birth and death rates (m in terms of the variable x become 

, {1 — w + w[ax + b{l — x)]} x{l — x) {1 — w + wlcx + d{l — x)]} x{l — x) 

-^^^^ " l-iv + w [ax^ + {b + c)x{l -x) + d{l - x)^] ' ^^^' " 1-w + w [ax^ + (b + c)x{l - x) + d(l - x)^] ^ ' 
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FIG. 1: (Color online). Shown is Invr,! as a function of n for a = 0.1, b = 0.7, c = 0.6, d — 0.2, w — 0.4 and A'^ = 150, so 
that n* = 75. The dynamics is implemented according to the IMP (M). We compare the analytical result (pOl) (solid line) with 
the numerical solution of the master equation (0) (dashed line) with rates (pq), and with a Gaussian approximation of this 
distribution with mean n* and standard deviation a = ^^ N / S" (x*) (dash-dotted line). Near the tails, one can clearly observe 
the non-Gaussian nature of the QSD. Note that, very close to the boundaries at n = 0(1) and N — n — 0(1) (see text), Eq. (Isd) 
has to be replaced by Eqs. (EOl) and (pi|). 



and one can check that 7^(0) = |7^(1)| = 1. The action S{x) is computed from Eq. (|T^) with the rates (|2^), yielding 



S{x) 



In 



1 — w + w[cq + d{\ — q)] 
1 — w + w[aq + 6(1 — q)] 



dq . 



(27) 



For further analytical treatment, it is convenient to introduce the following parameters (also used in Rcf. p3[): 
A = 1 — w + wa, B = 1 — w + wb, C = 1 — w + wc, and D = 1 — w + wd, where for ACG, C > A and B > D. 
Performing the integral (p7|), one obtains after some algebra 



e-^s(-) ^ [Ax + B{1 - x)]^^-^(^^) [Cx + D{1 - :,)]-^--^(c^) 

It can also be checked that 

BC-AD 



S"{x) = 



[Ax + B(l - x)][Cx + D(l - x)] 



is positive over the entire region < x < 1. It therefore follows from (|19|), and Eqs. (P6|), (P8|), (|2 
i.e. for A^-i/2 ^ X < 1 - N-^'"^, the QSD reads 



(28) 

(29) 
that in the bulk, 



7r(a;) = 



{C-A){B-D) 



[Ax + B{l~x)] 



N^-N(-B^) 



BC-AD 

^2ttN{BC-AD)T+{x)T-{x){B-A+C-D) [Cx + D{l-x)f''"'^'^'^^ \C-A+B-D 

Moreover, tti and tttv-i are obtained from Eqs. (p^): 

' (C - A){B - Df 



N{BC-AD) 
(B-A)(C-D) 



.(30) 



N 



TTl 



7rjv_i 



2t:BD{BC - AD) C-A + B-D \C -A + B - D 



BC-AD 



N{BC-AD) 
{B-A)(C-D) 



N 



{C-Af{B-D) 



2t:AC{BC-AD) C -A + B - D \C -A + B -D 



BC-AD 



N{BC-AD) 
{B-A){C-D) 



B-^{B^)D-^ic^) 
A-^(iAi)c-^(c^). (31) 



Eq. (|30|) determines the QSD for ACG evolving according to the IMP [close to the boundaries, one must use Eqs. ( pOJ ) 
and ( plD instead]. Clearly, the resulting QSD is non-Gaussian, which is especially evident near the tails. A typical 
example is shown in Fig. m where excellent agreement is observed between our analytical results and a numerical 
solution of the corresponding master equation (0). It is worth noticing that our theory is applicable when iV[S'(l) — 
S{x*)] ^ 1 and N[S{0) — S{x*)] ^ 1, which imposes a lower bound on w. Hence, while it is inapplicable in the 




FIG. 2: (Color online). Shown is Inr" as a function of the population size A'^, for a — 0.1, b = 0.7, c = 0.6, d = 0.2, with 
w = 0.2 in panel (a) and w = 0.7 in panel (b). Excellent agreement is observed between the analytical solution (solid line), 
given by Eqs. (psl) and (pll), and the numerical solution of the master equation (x's). 




FIG. 3: (Color online). Shown is lnr~^ as a function of the selection strength w, for a — 0.1, b — 0.7, c = 0.6, d — 0.2, with 
TV — 100 in panel (a) and A'^ = 200 in panel (b). Excellent agreement is observed between the analytical solution (solid line), 
given by Eqs. (bsl) and (pll), and the numerical solution of the master equation (x's). 



weak selection limit Nw ^ 1, recently investigated by other techniques (see e.g. [p3| , |30| and references therein), our 
approach successfully applies to the more general case of finite selection intensity < w < 1. 

The MFTs can now be found by using Eqs. (H) and (JH), with T1{0) = |T+(1)| = 1. As illustrated in Fig. |, 
the unconditional MFT asymptotically exhibits an exponential dependence on the population size N. That is, 
T ex iV^/^e^[^^"^^^ ^1, where S = min [5(0], S'(l)] depends on the entries of the payoff matrix and on the selection 
intensity w. It follows from Eqs. (25), (Eq) and ( |3l| ) that in the biologically relevant case of small (but not too 
small) selection intensity N^^ < w < 1, the MFTs grow exponentially as r^ ~ jsfi/2 ^Nw{a-c)y[2{c-a+b-d)]^ ^b ^ 
j^i/2^Nw{b-df/[2{c-a+b-d)]^ j^^j ^ ^ ^A^B |(^^AJ^^B^^ ^ ^^^^^A^^By j^ ^j^g opposite limit of large selection strength 
If — ^ 1, one can show that our results in the leading order coincide with those of Ref. [E4|. For finite selection strength, 
the exponential dependence of r is found to increase monotonically with w, as shown in Fig. y. Note, that in this 
figure and in all other figures (except Fig. ph, when simulating the master equation (ph, the initial number n of A's 
was chosen to be sufficiently large to avoid immediate fixation prior to reaching the metastable state. 

The ratio cf)^ /(f)^ — cj)^^ — \ allows to assess the infiuence of selection by comparing the fixation probability of 
species A and B for finite w. It follows from (11) that the fixation probability ratio (f> /cj)^ is given by 
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FIG. 4: (Color online). Shown is the ratio (f) /(j> of the fixation probabihties of species A and B as a function of the selection 
intensity w in the fMP. The theoretical prediction [Eq. (p2)] (solid) is compared with numerical solution ( x 's) of the master 
equation (pi). The parameters are a — 0.1, b = 0.7, c = 0.7 and d = 0.2, with A'' — 100 in panel (a) and A'' = 300 in (b). 



In Fig. M, the asymptotic expression (p2) as a function of the selection strength is compared with the numerical 
solution of the master equation (m, demonstrating an excellent agreement. This figure illustrates the exponential 
dependence of the ratio (j)^ /(p^ on w with a marked nonlinear behavior of the exponent. One can see a steep decay as 
the selection's strength increases (the fixation of B is thus more likely), and for w close to 1 the fixation of A is almost 
improbable. In the neutral case arising when w = 0, the stochastic dynamics is diffusive and the ratio of the fixation 
probabilities then strongly depends on the initial number n of A's. In stark contrast, for finite w cj)"^ /(p^ becomes 
independent of n (provided that n ^ 1 [H) and converges towards Eq. (^2|). This is demonstrated in Fig. 1^, where 
we have compared Eq. (53) with the numerical solution of the master equation (ph for various initial conditions. Note, 
that for w > and N ^- oo the fluxes to the absorbing states are vanishingly small and (^) becomes singular, with 
(j)^/(t>^ -> or (j)^/(p-^ — > oo depending on the rest of the parameters. When w — ^ 0, one has ip^/cj)^ — > x/{l — x). 



2. Linear Moran and local update processes 

The cases of the LMP and LUP, with rates T±{x) obtained respectively from Eqs. (P) and (ph, can be studied in 
the same manner as the fMP. Given the birth and death rates 7±(a;), one obtains the action [see Eq. (|6|)] and, with 
Eqs. (Ill), (p^), and (p^), one can calculate the QSD, fixation probability and MFTs. This leads essentially to the 
same qualitative features as in the fMP with low selection intensity w. Our findings are summarized in Fig. pi where 
our prediction for the unconditional MET for the LUP is found to grow exponentially with N and w, in excellent 
agreement with numerical results 



3. Fermi process 



We now consider the FP that has recently received considerable attention (see e.g. psl ESl pOl). As above, the 
transition rates 7± {x) for the FP are obtained from Eq. (^ . With Eq. ( p^ , the action in this case is quadratic 



S{x) 



and S"{x) = w{c - a + b- d) > 0. Using Eqs 
for TTi and ttn-i- 



w[{c - a)q + {d~ b){l - q)]dq = wx{d - b) [1 - x/{2x*)] , (33) 

j) and (p^) , after some algebra one obtains the following expressions 



TTl 



ttat-i = 



Nw {c-a){b-d) . ("-")"- 

ZTT [c— a + b — a)'V^ 



Nw (c—a)(b-d) . , ,, , , (a-a)^ 



exp 



wN{b - df 
2{c- a + b~d) 

wN{c-af 



2{c-a + b-d) 



(34) 
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It follows from these results that in the case of the FP, our approach is valid as long as w ^ A''"^, with N ^ 1. 
Therefore, our results are complementary to those of earlier works on this model which were carried out in the weak 
selection limit by treating selection as a linear perturbation to the neutral case w = (see, e.g., Refs. |23, E3, [30|). 
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FIG. 5: (Color online). Shown is the fixation probability of the A species versus the selection intensity w for the fMP. The 
dashed line is the analytical prediction given by Eqs. (U) and (p4|). The four solid lines are numerical solutions of the master 
equation (H) starting with n = 5, 10, 20, 30 (bottom to top) initial A's. The numerical results are found to converge towards 
the analytical prediction when w increases, with a convergence that improves when n increases. Inset: Ratios of the above four 
numerical curves and the analytical prediction (top to bottom). One clearly observes that the smaller w is, the larger n must 
be to avoid fixation prior to reaching the coexistence state. The parameters are a = 0.1, b = 0.7, c = 0.6, d = 0.2 and N = 300. 
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FIG. 6: (Color online), (a) For the LUP (m), shown is lnr~^ as a function of the population size A'^ and parameters w = 0.5, 
a = 0.2, h — c — 0.9, and d = 0.1. In (b), shown is Inr^^ for the LUP as a function of the selection intensity w, for the same 
process and same parameters as in (a) but A'^ = 100. An excellent agreement between our analytical predictions (solid lines) 
and the numerical solution of the master equation (x's) is observed. 

Using Eqs. (|l^), (p5|), ( p3| ) and (p4|) one can explicitly obtain the QSD, the fixation probability and the (uncondi- 
tional and conditional) MFTs. Here, one obtains the following asymptotic behavior of the conditional MFTs: 



iVl/2 



exp 



wN[c - af 



2(c-a 



while the unconditional MFT satisfies r 



T T 



b-d) 



Ari/2 



exp 



wN{b-dy 



2{c- a + b- d) 



(35) 



T^). We thus notice that, to leading order, the MFTs for the 



FP coincide with those obtained for the fMP in the limit of small (but not too small) selection strength A^ ^ ^ w <C 1. 
In addition to the MFTs, an interesting quantity to compute is the ratio of the fixation probabilities (j)^ and 0^ 

exp . (36) 



(p sinh[(c — a)w] 



(j)^ sinh[(6 — d)w 
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This ratio is larger than unity ii b — d > c — a, i.e. when x* > 1/2, which simply means that the closer x* is to 1, the 
easier it is for species A to fixate. In Fig. |7| we show the ratio between cj)-^ and (f)^ given by Eq. (^6|) as a function 
of N for a low and high selection strength {w — 0.2 and w — 0.7, respectively) and find excellent agreement with the 
numerical solution of the master equation (|^) . One can easily show that our asymptotic result ( p6[) coincides in the 
leading order with the exact result found from Eq. (H) (which takes a simple form in this case), where the difference 
in the pre- factor stems from self- interaction terms that were not excluded in our treatment [see Eq. (0)] |P3| . 



150 




150 



FIG. 7: (Color online). The ratio (p /cj>^ as a function of A*' for the FP (M. Comparison between the analytical result given by 
Eq. teq) (solid line) and the numerical solution of the master equation (ph (x's). The parameters in (a) are a — 0.5, b = c — 2, 
d = 0.4, i.e. X* ~ 0.516, and w = 0.2, and in (b) a = 0.3, h = c = 1, d = 0.2, i.e. x* ~ 0.533, and w — 0.7. As x* > 1/2, we 
notice that the fixation probability of species A increases with w and N. 



III. COORDINATION GAMES: FIXATION PROBABILITY 

In this section we use the WKB method to asymptotically compute the fixation probability in CG for N ^ 1 and 
finite w. After the presentation of the general treatment, our theoretical results are applied to the fMP (||) and FP (|^) 
update rules, and are compared with those obtained from the FPE. 



General treatment & results 



In CG, the elements of the payoff matrix satisfy d > b and a > c. Thus, within the realm of rate equations, the 
fixed point x* (g) is a repellor, whereas the absorbing states a; = and x — 1 are attractors. In the presence of noise, 
the fixation of the intruding species occurs rapidly p4| and therefore the main interest is in the fixation probability 
(/)„ : the probability that starting from n < n^, individuals of type A, the species A will fixate the population. In terms 
of the transition rates r='=(n), (j)^ obeys the following difference equation [0, p3, ^, p8| 



T+{r 



"n+l 



+ T-{n) 



bti 



[T+{n)+T-{n)]^^=0, 



(37) 



with boundary conditions (J)q ~ 0, 0^ = 1. Here, the probability (/>„ that the A's fixate starting from n individuals of 
type A is given by a sum of two components. The first is the probability to fixate starting from n -\- 1 A's multiplied 
by the probability to jump to state n + 1 from state n, which is T+(n)/[r+(ri) -I- T~{n)]. The second component is 
the probability to fixate starting from n ~ I A's multiplied by the probability to jump to state n — 1 from state n, 
which is r^(ri)/[T+(n) + T^{n)]. Note, that in this section (and differently from the treatment of ACG), as there is 
no metastability, the results strongly depend on the initial condition, that is, on the initial number of A's. 

As (j)„ = (j) (x) is a cumulative distribution function, it is more convenient to consider the corresponding PDF, 
defined as Vn = Vix) = (1)^+1 — <!>„■ 'Pi^) is peaked in the vicinity of x* = n^jN (see insets of Figs, g) and can be 
shown to satisfy Vq — 
r-(x)[0^(x)-0^(a;- 



^t 



Vn-1 = l-f^N-i andJ2o~^'Pn = 1. Rewriting Eq. (|3) asT+{x)[(t,'^{x + N~^)-<j)'^{x)]- 



N ^)] = and using the definition of 'P{x), one obtains the following equation for 'P{x) 
T+{x)V(x) - T-{x)V{x - N-^) = . 



(38) 
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This equation is similar in form to the QSME ( |l^ ) and is treated using the WKB ansatz 

V{x) ^ AcG e^^'5(^)-'Si(-). (39) 

Substituting (|39|) into (|38|), one obtains in the leading 0(1) order 7+(x) — T-{x)e^ ^^'> = 0, whose solution reads 



Six) = -S{x) 



HT+iO/T-iOH- 



(40) 



Here, S{x) is the negative of the expression (16), and thus S"{x*) < 0. In the sublcading 0{1/N) order, one obtains 
Si{x) = {l/2)S'{x) = (1/2) ln[7+(a;)/71(a;)]. As in Sec. II, the constant ^cG in Eq. (p9|) is found by normalizing the 
Gaussian asymptote of Vix) in the vicinity of x* . With Eq. (EQ), the final result for V{x) reads 



Vix) = 




,NlS{x)-S{xn]_ 



(41) 



As shown in the appendix, this result coincides up to subleading-order corrections with the exact solution of Eq. (|3S 
In particular, this implies that the recursion solution of Eq. ( p8[) near the boundaries a; = and x — 1 exactly coincides 
with the WKB result at a; ^ 1 and 1 — x <^ 1, respectively, and no special treatment is thus required near those 
boundaries (see also Ref. |3^). With Eq. ([4l[), we can write down the fixation probability, cj}^ — 'J27=o ^i' ^^ 




^NlS{i/N)-S{xn] 



(42) 



|), where its accuracy holds in 



This expression gives the fixation probability of species A for any finite w (see Figs, 
the entire region of x including the boundaries. 

Remarkably, the summation in Eq. (B2[) can be drastically simplified for x <^ x* , i.e., n ^ritf. This corresponds to 
the biologically important limit of the fixation probability of a few intruders of type A in a sea of B's |0, 0, IsO] . It is now 
convenient to split our discussion into two cases. For small (but not too small) selection intensity N^"^ ^ w ^ 1, P(a;) 
is slowly varying, and the sum (42 ) can be safely transformed into an integral. In this case, the term •\/7^(a;)/7+(x) ~ 1 
can be omitted from the integration, yielding 



0^(x) 



N\S"{x*)\ r N[S{ei-S{x')\ 



27r 



de 



(43) 



This approximation is valid when {Vn+i/Vn — 1| <C 1, i.e., when |S"(a::)| ^ 1, which assures that T+{x) ~ T~{x) and 
holds excellently for w ^ 1. 

In the second case, w = 0(1), the sum ( [l2|) is dominated by its last term when 1 ^ n <C n*. In this case, 
denoting k = n — 1 — i, one can Taylor-expand the summand about i = n in E q. ([4^ ), yielding e-^'S(i/7V)-5i(j/Af) ^ 
^-NS(n/N)-s^(n/N)+(k+i)s'(n/N)+Oii/N)^ Plugging this expression into the sum (|4|)7one has [with S'{x) = -S'{x)] 



n-l 



\x)^V{x)J2, 



~{k+l)S'ix) 



V{x) 



fe=0 



.S'(x) _ y 



(44) 



where S'{x) > for x < x*, and we have replaced the upper limit of the sum by infinity. Results (43) and (B3) are 



valid for N" <t^ x <t^ x*. Clearly, similar approximations can be made near a: = 1, in the region N^^^ 1 — a; <C 1. 

In Fig . B , for the fMP, we compare the numerical results for 0"^ with the WKB solution (E2) and with its approxi- 
mation (H4), and an excellent agreement is observed. The latter improves as w is increased. From (44), we infer that 
the fixation probability is exponentially small for finite selection intensity and therefore one generally has (j> {x) < x 
when X <$^ 1. In stark contrast with the weak selection limit [0, O, p3, ROl pM, this implies that when w is finite, 
selection always opposes the replacement of the wild species (B) by mutants (A). 



B. Applications 



We now apply the above general results to the cases of CG evolving according to the fMP ( 
compare our theoretical results with those of the FPE. For the fMP, the rates are given by {[ 



) and FP (|7D and 
with A > C and 
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FIG. 8: (Color online). Fixation probability of species A evolving according to the fMP m): comparison between the expression 
(tt2|) (dashed), its approximation (Uj) (x's), and numerical results (solid) as function of a;. Note that the theoretical (U2) and 
numerical results are indistinguishable. The parameters are A'' — 100, a — 1.2, 6 = 0.1, c = 0.3, d = 1.1, with w = 0.2 and 
w = 0.7 in the left and right panels, respectively. The quality and range of validity of the approximation (k4|) improves as w 
is increased. In the left panel the small-w approximation ( |43[ ) is also reported (dash-dotted). In this case, as S'{x) ^ 1, the 
approximation {nM is superior to (M), see text. In the insets of both panels, we show a comparison between the analytical 
[Eq. ([ll|)] (dashed) and numerical (solid) results for V{x), and excellent agreement is found over the entire range < a; < 1. 



D > B. In this case the action S{x) is given by Eq. ( p8| ) 
n = Nx <C n* individuals, is given by Eqs. (Esh and (M), i.e. 



The fixation probability of species A starting with 



0^(x) 



N 



VmpiOd^ 



"PfMp(a;) 



for TV ^ < w < 1 
for w = 0(1), 



(45) 



where 7^fMp(a^) is given by Eq 
the theoretical predictions [Eq. 
< X < 1 . Results (^5|) generalize the results of Ref . 
and by including the subleading-order correction. 
For the FP (R) starting with n ^ n* mutants, Eqs 



l|) and S{x) given by (|28|). In Fig. |g, we compare (for w — 0.2 and w — 0.7) 

)] with the numerical results and find an excellent agreement over the entire range 

by considering arbitrary (finite) selection strength < w < 1 



cf>^{x) 



3|) and (^4|) can be explicitly calculated. Using (^3|) one finds 
- jerf y/Nx*a{x/x* - 1) + eri{y/Nx*a)\ for iV"^ < w < 1 

for w = 0(1). (46) 



a 

nNx* 



-a{x/x'-l)[N{x~x*) + l] 



^~2a{x I x' -\) _ Y 



where erf(x) = (Ij^) Q e'"^ dy denotes the usual error function, and a — w{d — h)/2 > 0. The small-w result 
coincides with the exact results in the continuum limit (see, e.g., Refs. ||2^, |30|). 

Fixation probabilities are often computed using diffusion approximations, like the FPE pBl p9|-^, that can be 
obtained from a van Kampen size expansion. This expansion implicitly assumes that </) (xjvarics slowly over the 
entire range of < x < 1 [E^ p8| . Here, we are interested in comparing the predictions of the FPE with those of the 
WKB approach when the selection intensity is small (but not too small) N~^ <C ui ^ 1. In this case, the fixation 
probability in both the WKB and Fokker-Planck treatments ||2^, ^, |3^ , is conveniently rewritten as 



0^(a:) = ^iM where *(a;) = / e^ lo ^i-)d^dC 
*(1) Jq 



Here, by comparing Eq. (pTJ) to Eq. 



^), one has 6wkb(2;) = A^ln [7+(a;)/71(a;)], while 

-T+ix)-T-ixy 



OfMx) - 2N 



T+{x)+T-{x) 



(47) 



(48) 
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FIG. 9: (Color online). Fixation probability (x) evolving according to the flVIP M) versus w. comparison between the 
WKB result [Eq. (||)] (solid), results of the full [Eqs. @ and @] (dashed) and linearized [Eqs. ^) and (^)] (dash-dotted) 
FPE, and numerical solution of Eq. (jsl) (x's). Parameters are a — 2, b = 0.2, c — 0.3, d — 1.8, A'' = 150, with initial 
condition x — n/N = 10/150. The WKB solution agrees excellently with the numerical results, while the FPE approximations 
systematically deviate from the WKB result as w increases. Inset: ratios between the WKB and the results of the full (dashed), 
and linearized (dash-dotted) FPE as function of w. The linearized FPE has a narrower region of applicability, see text. 



Furthermore, the FPE is often considered within the linear noise approximation [E7l 
(p8|) is expanded to linear order in a; — x* , yielding 



e 



£FPE 



(x) = 2N{x - X*) 



T[{,x*)-TL{x*) 



T+{x*) + T-{x* 



M. In this case, 8fpe(2:) 



(49) 



To see how the WKB approximation compares with that of the FPE, we Taylor expand the functions 6 around a;*, 
and compute A9FPE(a^) — 0WKB(a^) — 0FPE(a;) and A0fFPE(a;) = ©wkbI^;) — QiYVEix). For the IMP one finds 



*\2 



AOfpeCx) ~ Cfpe Nw%x -x*Y , ABfFPEla;) ~ C^fpe Nw\x - x*) 



(50) 



where Cfpe = {lll2){a - h - c+ df /[{a - h - c + d){l - w) + {ad - hc)w^ and Qfpe ^ il/2){a-b + c- d){a -b - 
c + d)'^/[(a — b — c + d){l — w) + {ad — bc)w]'^ , are both 0(1). Results ( pO| ) demonstrate that the exponent 6 of our 
theory and those obtained from the FPE significantly deviate from each other when x — x* = 0(1) [e.g. when x -^ 1 
or 1 — X <C 1] and w is finite. Looking at Eq. ( pG), the results of the full and linear FPE agree with the WKB theory 
and numerical calculations (see also Figs. |9| and |10|) when N^^ <t^ w <$^ N^^^^ and N^^ ^ w ^ N~^^^, respectively. 
This implies that demanding that w <C 1 does not guarantee the applicability of the FPE, as the results of the full 
and linear FPE are plagued by exponentially large errors already when w > N^^^^ and w > N'^^"^, respectively. 
While the predictions of the FPE further deteriorate when w increases, our theory improves and allows the accurate 
calculation of the exponentially small fixation probability for any finite value of w (see Fig. p|)p5[. 



IV. SUMMARY AND CONCLUSION 



In this work, we have studied large- fluctuation- induced fixation in (2x2) anti-coordination and coordination evolu- 
tionary games using a WKB-based approach. In both classes of games, the deterministic description is characterized 
by the existence of an interior fixed point separating two absorbing fixed points. The latter are the only possible 
outcomes of the dynamics when internal stochasticity is taken into account. Yet, for ACG, the mean time necessary 
to reach one of the absorbing fixed points (mean fixation time) is typically very large as the system typically spends 
an exponentially long time in the metastable coexistence state. On the other hand, in CG fixation occurs rapidly, 
and the quantity of interest is the (small) probability that a system comprising of a few mutants takes over the entire 
population, causing the extinction of the wild species. As the stochastic dynamics of evolutionary games is formulated 
in terms of one-dimensional single-step birth-death processes (with frequency-dependent rates), there exist exact for- 
mulas for the mean time and probability of fixation. However, these unwieldy expressions are nontrivial to analyze 
and cannot be generalized to multi-step processes or higher dimensions. To circumvent this difliculty, one popular 
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FIG. 10: (Color online). Shown are the ratios of the WKB result for (j) (x) [Eg . (|45[ )] with the numerical result (solid), with 
the fuU FPE [Eqs^ @ and @] (dashed), and with the linearized FPE [Eqs. (|^and (||)] (dash-dotted), as a function of 
iV, for the fMP @. Here x = 10/iV (i.e. n = 10), a = 2, 6 = 0.2, c = 0.3, d = 1.8 and w = 0.5. The WKB predictions agree 
excellently with the numerical results, while there are systematic deviations (that increase with A^) between the latter and the 
predictions of the the diffusion approximations (full and linearized FPE), see text. 



approach is to resort to diffusion approximations, as those based on the FPE. However, the latter are ill-suited to 
describe phenomena like fixation that result from large fiuctuations. 

Here we have presented an alternative approach based on the WKB theory, which we have generalized to account 
for the existence of multiple absorbing states. Within this approach, the stochastic dynamics of the system is formally 
mapped, in the leading order, onto a Hamiltonian system, whose nontrivial zero-energy trajectory encodes, with 
the maximum probability, the rare event in question. By using the WKB approach complemented by recursive 
solutions of the master equation near the two absorbing boundaries, we have obtained general results, for the complete 
statistics including large fluctuations, of generic one-dimensional birth-death systems possessing two absorbing states. 
Our results were obtained including important pre-exponential factors, which were found to scale as some power of 
the typical population size. Along with the generic treatment, we have also considered the most frequently used 
microscopic dynamics, based on the frequency-dependent Moran process (fMP), the linear Moran and local update 
processes, as well as the Fermi process (FP). In particular we have focused on the fMP and FP and have obtained 
explicit analytical results for the complete metastable probability distribution function of population sizes and for the 
fixation probability and times in ACG, as well as for the fixation probability in CG. All these results were obtained 
for arbitrary and finite w, which allowed us to shed further light on the combined influence of selection and random 
fluctuations in evolutionary processes. Finally, by comparing our analytical results to those of the FPE, we have 
explicitly found the region of applicability of the FPE and have shown that the WKB approach is vastly superior over 
the FPE when the selection strength is finite. 

While we have here focused on 2 x 2 evolutionary games, the WKB-based method presented in this work can 
be generalized and is expected to be useful for multi-species problems, like the rock-paper-scissors games that have 
recently received considerable attention (see, e.g., m Esj). 



Appendix 

In this appendix, we show that the WKB result (Bl|) asymptotically coincides with the exact solution of (pS 
by using recursion, the exact solution of (pq) reads 



-pc 



= r. 



n 






Poexp 



±iJl^ 



yr+iz] 



Indeed 



(Al) 



Here, Vo is a constant to be found by normalization. For A^ ^ 1, the sum in the exponent of Eq. (Al) can be 
transformed into an integral using Euler-Maclaurin formula X]r=o /(*) ~ In /(2^)'^^+(l/2)[/(")+/(0)] + (l/12)[/'(n)-f 
/' (0)] — As the sum is in the exponent of ( Al ) , such a transformation should be done carefully and subleading-order 



16 



corrections to the integral must be taken into account. Therefore, for A^ ::^ 1 and x = n/N, one has 



Ein 



^^-«^.,,/^/^^,, la^i^), (A2) 



,^1 \T+{^)J- Jo \T+{oj^'2 Vrz(o)/r;(o);' 



up to 0{1/N) corrections, where we used the fact that T-iO)/T+iO) ^ T_^(0)/7I(0). Thus, Eq. (Al) becomes 



V {x)-Vo^^,^^^^^^^^^e . (A3) 

where we have used the definition of S{x) from Eq. (|40|). Finally, Vo is determined by demanding that ^q Vn = 1- 
As before, for N ^ 1 we approximate this sum by an integral, whose main contribution arises from the Gaussian 
region around x c^ x* , where the function V{x) varies slowly. By doing so, one obtains 



^" V 2TrN \Tm ■ ^ ' 

With this result, (A3) coincides with Eq. ( pl| ) obtained directly from our WKB treatment. 
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